%% 模拟退火解决不回到起点的TSP问题(很多地方我直接引用的蒙特卡罗模拟里面的代码),我找到的最佳是6615.3
tic
clear;clc
% 38个城市，TSP数据集网站(http://www.tsp.gatech.edu/world/djtour.html) 
coord = [11003.611100,42102.500000;11108.611100,42373.888900;11133.333300,42885.833300;11155.833300,42712.500000;11183.333300,42933.333300;11297.500000,42853.333300;11310.277800,42929.444400;11416.666700,42983.333300;11423.888900,43000.277800;11438.333300,42057.222200;11461.111100,43252.777800;11485.555600,43187.222200;11503.055600,42855.277800;11511.388900,42106.388900;11522.222200,42841.944400;11569.444400,43136.666700;11583.333300,43150.000000;11595.000000,43148.055600;11600.000000,43150.000000;11690.555600,42686.666700;11715.833300,41836.111100;11751.111100,42814.444400;11770.277800,42651.944400;11785.277800,42884.444400;11822.777800,42673.611100;11846.944400,42660.555600;11963.055600,43290.555600;11973.055600,43026.111100;12058.333300,42195.555600;12149.444400,42477.500000;12286.944400,43355.555600;12300.000000,42433.333300;12355.833300,43156.388900;12363.333300,43189.166700;12372.777800,42711.388900;12386.666700,43334.722200;12421.666700,42895.555600;12645.000000,42973.333300];
n = size(coord,1);  % 城市的数目
coord = [coord (1:n)']; % 新加一列表示城市的编号(等会儿打印路径有用)
% 设置起点城市和终点城市
begin_city = 5;  % 起点城市
end_city = 20; % 终点城市
tem = coord;  % 备份下原来的坐标
coord([begin_city end_city],:) = []; % 将起点和终点坐标先删除
coord = [tem(begin_city,:); coord; tem(end_city,:)];  % 重新调整坐标，此时起点在最开始，终点在最后面

figure  % 新建一个图形窗口
plot(coord(:,1),coord(:,2),'o');   % 画出城市的分布散点图
hold on % 等一下要接着在这个图形上画图的


d = zeros(n);   % 初始化两个城市的距离矩阵全为0
for i = 2:n  
    for j = 1:i  
        coord_i = coord(i,:);   x_i = coord_i(1);     y_i = coord_i(2);  % 城市i的横坐标为x_i，纵坐标为y_i
        coord_j = coord(j,:);   x_j = coord_j(1);     y_j = coord_j(2);  % 城市j的横坐标为x_j，纵坐标为y_j
        d(i,j) = sqrt((x_i-x_j)^2 + (y_i-y_j)^2);   % 计算城市i和j的距离
    end
end
d = d+d';   % 生成距离矩阵的对称的一面


%% 参数初始化
T0 = 1000;   % 初始温度
T = T0; % 迭代中温度会发生改变，第一次迭代时温度就是T0
maxgen = 200;  % 最大迭代次数
Lk = 500;  % 每个温度下的迭代次数
alfa = 0.95;  % 温度衰减系数

%%  随机生成一个初始解
path0 = randperm(n-2)+1;  % 生成一个2至n-1的随机打乱的序列
path0 = [1, path0, n]; % 生成初始的路径
result0 = calculate_tsp_d_homework(path0,d); % 调用我们自己写的calculate_tsp_d_homework函数计算当前路径的距离

%% 定义一些保存中间过程的量，方便输出结果和画图
min_result = result0;     % 初始化找到的最佳的解对应的距离为result0
RESULT = zeros(maxgen,1); % 记录每一次外层循环结束后找到的min_result (方便画图）

%% 模拟退火过程
for iter = 1 : maxgen  % 外循环, 我这里采用的是指定最大迭代次数
    for i = 1 : Lk  %  内循环，在每个温度下开始迭代
        path1 = gen_new_path_homework(path0);  % 调用我们自己写的gen_new_path_homework函数生成新的路径
        result1 = calculate_tsp_d_homework(path1,d); % 计算新路径的距离
        if result1 < result0    % 如果新路径的距离小于当前路径的距离
            path0 = path1; % 更新当前路径为新路径
            result0 = result1;
        else
            p = exp(-(result1 - result0)/T); % 根据Metropolis准则计算一个概率
            if rand(1) < p   % 生成一个随机数和这个概率比较，如果该随机数小于这个概率
                path0 = path1;  % 更新当前路径为新路径
                result0 = result1;
            end
        end
         % 判断是否要更新找到的最佳的解
        if result0 < min_result  % 如果当前解更好，则对其进行更新
            min_result = result0;  % 更新最小的距离
            best_path = path0;  % 更新找到的最短路径
        end
    end
    RESULT(iter) = min_result; % 保存本轮外循环结束后找到的最小距离
    T = alfa*T;   % 温度下降       
end

% 打印我们的路径，此时coord第三列就可以帮我们还原真实的编号了
path_str = ''; % 生成一个空的字符串
for i = 1:n-1 
    path_str =  strcat(path_str, num2str(coord(best_path(i),3)), '-->');  % 不断更新这个字符串(strcat是字符串拼接函数)
end
path_str =strcat(path_str, num2str(coord(n,3))); % 别忘了加上终点
disp('最佳的方案是：'); disp(path_str)
disp('此时最优值是：'); disp(min_result)


for i = 1:n-1 
     j = i+1;
    coord_i = coord(best_path(i),:);   x_i = coord_i(1);     y_i = coord_i(2); 
    coord_j = coord(best_path(j),:);   x_j = coord_j(1);     y_j = coord_j(2);
    plot([x_i,x_j],[y_i,y_j],'-b')    % 每两个点就作出一条线段，直到所有的城市都走完
%     pause(0.02)  % 暂停0.02s再画下一条线段
    hold on
end

%% 画出每次迭代后找到的最短路径的图形
figure
plot(1:maxgen,RESULT,'b-');
xlabel('迭代次数');
ylabel('最短路径');

toc
